function OOSStats = OOS(X,Y,Rf,MktVar,h,Tbeg,Tend,OOSbeg,gam)

load ClarkMcCrackenCriticalValues;
T = length(Y);
R = OOSbeg - Tbeg - h; P = Tend - OOSbeg; pi = P/R;
K = 1;

YMdl = NaN(T,1); YUnc = NaN(T,1);
WtMdl = NaN(T,1); WtUnc = NaN(T,1);
for t = OOSbeg-1:Tend-h 
    % the first forecast is for OOSbeg. 
    % If h=3 (quarerly) then the first forecast will be for 194603 (if OOSbeg=194601)
    % If h=12 (annual) then the first forecast will be for 194612 (if OOSbeg=194601)  
    x = X(Tbeg:t-h,t); y = Y(Tbeg+h:t);
    x = [ones(length(y),1) x];
    bols = x\y;
    YMdl(t+h) = [1 X(t,t)]*bols;
    YUnc(t+h) = mean(y);
    
    WtMdl(t+h) = YMdl(t+h)/(gam * MktVar(t));
    WtUnc(t+h) = YUnc(t+h)/(gam * MktVar(t));
end
ErrMdl = Y - YMdl; ErrUnc = Y - YUnc;
EMdl = ErrMdl(OOSbeg+h-1:Tend); EUnc = ErrUnc(OOSbeg+h-1:Tend);
OOSStats.ErrMdl = ErrMdl; OOSStats.ErrUnc = ErrUnc;

SE1 = EUnc.^2; SE2 = EMdl.^2;
d = SE1 - SE2;
c = EUnc.*(EUnc - EMdl);
OOSStats.R2 = 1 - mean(SE2)/mean(SE1);
OOSStats.DMAE = mean(abs(EMdl)) - mean(abs(EUnc));
OOSStats.MSET = [sqrt(P)*mean(d)/std(d) interp1(Crtpi,MSE_T_90(K,:),pi) interp1(Crtpi,MSE_T_95(K,:),pi) interp1(Crtpi,MSE_T_99(K,:),pi)];
OOSStats.MSEF = [P*mean(d)/mean(SE2) interp1(Crtpi,MSE_F_90(K,:),pi) interp1(Crtpi,MSE_F_95(K,:),pi) interp1(Crtpi,MSE_F_99(K,:),pi)];
OOSStats.ENCT = [sqrt(P)*mean(c)/std(c) interp1(Crtpi,ENC_T_90(K,:),pi) interp1(Crtpi,ENC_T_95(K,:),pi) interp1(Crtpi,ENC_T_99(K,:),pi)];
OOSStats.ENCF = [P*mean(c)/mean(SE2) interp1(Crtpi,ENC_F_90(K,:),pi) interp1(Crtpi,ENC_F_95(K,:),pi) interp1(Crtpi,ENC_F_99(K,:),pi)];
OOSStats.ENClambd = (EUnc-EMdl)\EUnc;

% Add Campbell/Thomson restriction - if forecast < 0, set it to zero
YMdl(YMdl<0) = 0; YUnc(YUnc<0) = 0;
ErrMdl = Y - YMdl; ErrUnc = Y - YUnc;
EMdl = ErrMdl(OOSbeg+h-1:Tend); EUnc = ErrUnc(OOSbeg+h-1:Tend);
OOSStats.ErrMdl_CT = ErrMdl; OOSStats.ErrUnc_CT = ErrUnc;

SE1 = EUnc.^2; SE2 = EMdl.^2;
d = SE1 - SE2;
c = EUnc.*(EUnc - EMdl);
OOSStats.R2_CT = 1 - mean(SE2)/mean(SE1);
OOSStats.DMAE_CT = mean(abs(EMdl)) - mean(abs(EUnc));
OOSStats.MSET_CT = [sqrt(P)*mean(d)/std(d) interp1(Crtpi,MSE_T_90(K,:),pi) interp1(Crtpi,MSE_T_95(K,:),pi) interp1(Crtpi,MSE_T_99(K,:),pi)];
OOSStats.MSEF_CT = [P*mean(d)/mean(SE2) interp1(Crtpi,MSE_F_90(K,:),pi) interp1(Crtpi,MSE_F_95(K,:),pi) interp1(Crtpi,MSE_F_99(K,:),pi)];
OOSStats.ENCT_CT = [sqrt(P)*mean(c)/std(c) interp1(Crtpi,ENC_T_90(K,:),pi) interp1(Crtpi,ENC_T_95(K,:),pi) interp1(Crtpi,ENC_T_99(K,:),pi)];
OOSStats.ENCF_CT = [P*mean(c)/mean(SE2) interp1(Crtpi,ENC_F_90(K,:),pi) interp1(Crtpi,ENC_F_95(K,:),pi) interp1(Crtpi,ENC_F_99(K,:),pi)];
OOSStats.ENClambd_CT = (EUnc-EMdl)\EUnc;

% Do utiity calculations. Weights are constrained to be between 0% and 150%
% These calculations use NON CT forecasts
WtMdl(WtMdl<0) = 0; WtMdl(WtMdl>1.5) = 1.5; RportMdl = WtMdl.*Y + Rf; RportMdl = RportMdl(OOSbeg+h-1:Tend);
WtUnc(WtUnc<0) = 0; WtUnc(WtUnc>1.5) = 1.5; RportUnc = WtUnc.*Y + Rf; RportUnc = RportUnc(OOSbeg+h-1:Tend);

UtilMdl = mean(RportMdl) - 0.5*gam*var(RportMdl);
UtilUnc = mean(RportUnc) - 0.5*gam*var(RportUnc);

V = cov(RportMdl,RportUnc);
Sig = [V zeros(2); zeros(2) 2*(V.^2)];
D = [1; -1; -gam/2; gam/2];
tDiffUtil(1) = (UtilMdl - UtilUnc)/sqrt(D'*Sig*D/P);

% add always 100% in market
RportMkt = Y + Rf; RportMkt = RportMkt(OOSbeg+h-1:Tend);
UtilMkt = mean(RportMkt) - 0.5*gam*var(RportMkt);
V = cov(RportMdl,RportMkt);
Sig = [V zeros(2); zeros(2) 2*(V.^2)];
D = [1; -1; -gam/2; gam/2];
tDiffUtil(2) = (UtilMdl - UtilMkt)/sqrt(D'*Sig*D/P);

% add 50% in market
WtMktG = 0.5;
RportMktG = WtMktG*Y + Rf; RportMktG = RportMktG(OOSbeg+h-1:Tend);
UtilMktG = mean(RportMktG) - 0.5*gam*var(RportMktG);
V = cov(RportMdl,RportMktG);
Sig = [V zeros(2); zeros(2) 2*(V.^2)];
D = [1; -1; -gam/2; gam/2];
tDiffUtil(3) = (UtilMdl - UtilMktG)/sqrt(D'*Sig*D/P);

OOSStats.WtMdl = WtMdl; OOSStats.WtUnc = WtUnc;
OOSStats.UtilMdl = UtilMdl; OOSStats.UtilUnc = UtilUnc; OOSStats.UtilMkt = UtilMkt; OOSStats.UtilMktG = UtilMktG; 
OOSStats.tDiffUtil = tDiffUtil;
return